% ME452 Project 2 bearing selection
% design must include 6 total bearings:
%   2 bearings on input shaft
%   2 bearings on output shaft
%   2 bearings on intermediate shaft

clear,clc

% DESIGN REQUIREMENTS
Rd = 0.99; % design reliability
Ld = 10; % desired life (continuous operation) [years]
n_f = 1; % minimum bearing fatigue factor of safety
n_s = 1; % static crush factor of safety, C0/(overload condition (3 x typical load))

% INPUT SHAFT CONDITIONS
n_A = 200; % input shaft speed [rpm] typical load
F_l = 260; % lateral load, weight of blades and hub [lb] typical load
F_a = 50.25; % axial load, from thrust [lb] typical load

% FBD RESULTS
% forces lbs, n_A rpm
switch '2'
    case '1' % BEARING 1 %-((-260*(1.05+1.05+3.75))+(82.58*(1.05)))/(1.05+1.05),-(683+82.5-260)
n_A = 200; F_r = 683; F_a = 0; F_r = abs((82.58*1.4455+-260*(1.4455+2.4695+4.61))/(1.4455+2.4695)),
    case '2' % BEARING 2 - supports axial load % 506
n_A = 200; F_r = abs(-(535.6655+82.58-260)), F_a = 50.25;
    case '3' % BEARING 3
n_A = 200*49/17; F_r = 70; F_r = abs(-(-24.75*(0.8145)-82.6*(.8145+4.105))/(.8145+4.105+1.334)), F_a = 0; % -(-24.75*(0.8)-82.6*(.8+3.2))/(.8+3.2+1.05)
    case '4' % BEARING 4
n_A = 200*49/17; F_r = abs(-(68.2033-24.8-82.6)), F_a = 0; %-(70-24.8-82.6)
    case '5'% BEARING 5
n_A = 200*49/17*104/17; F_r = 24.75/2, F_a = 0; 
% BEARING 6 % identical loading to 5, symmetry
    case '6'
        % test values from class
n_A = 60e6/(Ld*365.25*24*60); F_r = 1000; F_a = 600; Rd = .995;  % switch af to 1.0 from 1.2
end

Ld_r = Ld*365.25*24*60*n_A; % desired life in revolutions

% first (input) bearing must handle the axial load
% iterative selection procedure, first solve for C10 with no axial load
a1 = 0.02 + 4.439*(log(1/Rd))^(1/1.483); % weibull parameters
af = 1.2; % commercial gearing p.583
C10req = af*F_r*(((Ld_r)/1e6)/(a1))^(1/3); % p. 578

BearingFe_table = load('BearingFE.mat','-ascii'); % load bearing table
for iter = 1:5 % should not take more than 5 iterations
    C10 = input(['C10 of smallest bearing with C10 greater than: ',num2str(C10req),' lbf:  ']);
    
    % check if axial load can be neglected, could make it simpler and
    % interpolate..
    ind = find(BearingFe_table(:,1)<=max(F_a/C10,.014));
    e = BearingFe_table(ind(end),1);
    if F_a/(1*F_r) <= e
        disp('neglect axial load')
        F_e = F_r;
    else
        disp('using equivalent loading')
        XX = BearingFe_table(ind(end),3);
        YY = BearingFe_table(ind(end),4);
        F_e = F_a*YY + 1*F_r*XX; 
        F_e
    end
    
    % solve for new C10
    C10req = af*F_e*(((Ld*365.25*24*60*n_A)/1e6)/(a1))^(1/3); % p. 578
    
    if ~(C10req>=C10)
        disp(['use bearing with C10 = ',num2str(C10),' acceptable'])
        break % solution found, exit loop
    end
end

% check static loading
C0 = input(['C0 of selected bearing to check static loading lbf:  ']);
Fe = sqrt((3*F_r)^2 + (3*F_a)^2);

for iters = 1:5
    if Fe<=1*C0
        disp('bearing acceptable')
        break % solution found, exit loop
    else
        C0 = input(['C0 of next larger bearing to check static loading lbf:  ']);
    end
end

ACTUAL_LIFE_YEARS = ((C10/(af*F_e))^3*a1)*1e6/(365.25*24*60*n_A) 
ACTUAL_LIFE_HOURS = ACTUAL_LIFE_YEARS*365.25*24
FS_YIELD = C0/Fe